{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "973A8B4F35C443BAB0DC10F6FDE4C949",
        "id": "DDB780EA379C403A9F118EE52D12AD83",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "# <center> Lecture 12 : Poisson & Negative Binomial Regression </center>  \n",
        " \n",
        "## <center> Instructor: Dr. Hu Chuan-Peng </center> "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "D3FB5D67B79D4F8085299A04C606457B",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 回顾：  \n",
        "- Lecture 09: 简单线性模型 (在PyMC中建立模型、先验预测检验、后验评估)  \n",
        "- Lecture 10: 模型评估 (模型评估的原则、模型拟合、后验预测的准确性)  \n",
        "- Lecture 11: 模型比较 (多元线性模型的建立、模型比较的原则、PyMC中模型比较)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AEAC337E349248CBA4DDE6DB8A77B90B",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "![Image Name](https://cdn.kesci.com/upload/image/rkz1ehen1l.png?imageView2/0/w/960/h/960)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "77C2CCE8D31247159738B0A489903BCF",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**大作业的选题情况？**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "548BBEC803A3438F94818B37E3290969",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 引入  \n",
        "\n",
        "* 真实世界中，我们感兴趣的变量（因变量）并不仅包括连续数据(continuous data)，也可能是二分数据（如被试是否做出正确回答）或是计数数据(count data，回答正确的个数)。  \n",
        "\n",
        "* 目前，我们先前的课程中关注了基于正态分布的回归模型，但它并非适用于所有数据类型。  \n",
        "\n",
        "* 在本节课中，我们关注如何对**计数数据**进行分析。  \n",
        "\n",
        "* 同时，我们将展示，面对真实的数据，我们如何探索不同的模型从而找到相对更优的模型  \n",
        "\t* 线性回归(正态分布)  \n",
        "\t* 广义线性模型(泊松分布)  \n",
        "\t* 广义线性模型(负二项分布)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BA0D87559D6C496DA1D0DAE8AD01BEA3",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 本课内容  \n",
        "\n",
        "* 如何对**计数数据**进行分析。  \n",
        "\n",
        "* 如何探索不同的模型，从而找到相对更优的模型。  \n",
        "\t* 线性回归(正态分布)  \n",
        "\t* 广义线性模型(泊松分布)  \n",
        "\t* 广义线性模型(负二项分布)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "D3321EB2BB0843BE928EB32FD6F4E7F8",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**数据集介绍**  \n",
        "\n",
        "* 使用的数据集仍然来自来自Human Penguin Project (Hu et al., 2018, doi: 10.17605/OSF.IO/H52D3)，本节课中我们关注**压力**和**是否恋爱**对**社交网络大小**的影响。  \n",
        "\n",
        "* 在这份数据中，社交网络大小的定义为，在每两周内，个体的社交对象数量(社交对象包含父母、子女、同事等等)  \n",
        "\n",
        "\n",
        "> * 数据来源: Hu, C.-P. et al. (2018). Raw data from the Human Penguin Project. Open Science Framework. https://doi.org/10.17605/OSF.IO/H52D3  \n",
        "> * 压力量表来源：Cohen, S., Kamarck, T. & Mermelstein, R. A global measure of perceived stress. J. Health. Soc. Behav. 24, 385–396 (1983).  \n",
        "> * 社交网络测量：Cohen, S. et al. Social ties and susceptibility to the common cold. JAMA 277, 1940–1944 (1997)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "A650A9D814B540B4A343246855431DD7",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**🤔研究问题：如何通过模型形式化这些问题？**  \n",
        "\n",
        "**问题0:  哪种分布可以更好的描述社交网络变量（本节课主要内容）**   \n",
        "\n",
        "问题1：个体压力与社交网络大小的关系  \n",
        "\n",
        "> * 个体在面对压力时，通常会寻求他人的支持和安慰，以缓解压力的影响。因此，压力可能会扩大个体社交网络  \n",
        "> * 然而，个体在压力大时也可能会选择封闭自己，减小自己的社交网络，这也被认为是抑郁症产生的原因。  \n",
        "\n",
        "问题2：恋爱情况与社交网络大小的关系  \n",
        "\n",
        "> * 恋爱关系可能使个体与伴侣的社交圈子相互连接，进一步扩大其社交网络。  \n",
        "> * 相反，恋爱双方可能缩小自己的社交网络，以此避免影响彼此关系的风险因素出现。  \n",
        "\n",
        "**问题3：个体压力和恋爱情况与社交网络大小的关系（本节课主要内容）**  \n",
        "\n",
        "问题4：个体压力和恋爱情况对社交网络大小的交互作用  \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AAA0D256B433407B80A375FB55E92009",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**原始数据可视化**  \n",
        "\n",
        "* 社交网络大小的分布情况  \n",
        "\n",
        "* 压力与社交网络大小的散点图情况  \n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/s4o8026lgp.png?imageView2/0/w/960/h/960)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "_id": "9ECC4BF549D0407D979D03464EDA284F",
        "collapsed": false,
        "id": "A28CA4BA746A4C3B92C422D55A69634F",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "skip"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 导入 pymc 模型包，和 arviz 等分析工具 \n",
        "import pymc as pm\n",
        "import arviz as az\n",
        "import seaborn as sns\n",
        "import scipy.stats as st\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import xarray as xr\n",
        "import pandas as pd\n",
        "\n",
        "# 忽略不必要的警告\n",
        "import warnings\n",
        "warnings.filterwarnings(\"ignore\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "collapsed": false,
        "id": "05D0AECDA324489DA662A7F773DC73D1",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "skip"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 通过 pd.read_csv 加载数据 Data_Sum_HPP_Multi_Site_Share.csv\n",
        "try:\n",
        "  df_raw = pd.read_csv('/home/mw/input/bayes20238001/Data_Sum_HPP_Multi_Site_Share.csv')\n",
        "except:\n",
        "  df_raw = pd.read_csv('data/Data_Sum_HPP_Multi_Site_Share.csv')\n",
        "\n",
        "df = df_raw[df_raw[\"Site\"] == \"Tsinghua\"]\n",
        "\n",
        "#1 = \"yes\"; 2 = \"no\"\n",
        "df[\"romantic\"] =  np.where(df['romantic'] == 2, 0, 1)\n",
        "df[\"romantic_recode\"] =  np.where(df['romantic'] == 1, \"yes\", \"no\")\n",
        "\n",
        "#设置索引\n",
        "df[\"index\"] = range(len(df))\n",
        "df = df.set_index(\"index\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "D92F0D56217E4F338529B464E73C4A26",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### Linear regression model  \n",
        "\n",
        "首先，我们使用上节课学习过的正态回归模型来进行建模分析，虽然我们知道计数数据可能不适合正态分布，但是，正态分布的数据将能提供一些直觉，告诉我们，模型错在哪里。  \n",
        "\n",
        "$$  \n",
        " \n",
        "\\begin{array}{lcrl}  \n",
        "\n",
        "\\text{data:} & \\hspace{.05in} & Y_i|\\beta_0,\\beta_1,\\beta_2,\\sigma & \\stackrel{ind}{\\sim} N(\\mu_i, \\; \\sigma^2)  \\;\\; \\text{ with } \\;\\; \\mu_i = \\beta_0 + \\beta_1X_{i1} + \\beta_2X_{i2} \\\\  \n",
        "\\text{priors:} & & \\beta_{0}  & \\sim N\\left(20, 10^2 \\right)  \\\\  \n",
        "                    & & \\beta_1  & \\sim N\\left(0, 10^2 \\right) \\\\  \n",
        "                    & & \\beta_2  & \\sim N\\left(0, 10^2 \\right) \\\\  \n",
        "                    & & \\sigma   & \\sim \\text{Exp}(0.6)  .\\\\  \n",
        "\\end{array}  \n",
        "\n",
        "$$  \n",
        "\n",
        "* 我们快速运行模型定义、MCMC采样和后验预测分布的结果"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "collapsed": false,
        "id": "AD9EE50E737F4D409D44C99431A54A8C",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "with pm.Model(coords = {\"obs_id\": df.index}) as normal_model:\n",
        "\n",
        "    beta_0 = pm.Normal(\"beta_0\", mu=20, sigma=10)          #定义beta_0          \n",
        "    beta_1 = pm.Normal(\"beta_1\", mu=0, sigma=10)           #定义beta_1\n",
        "    beta_2 = pm.Normal(\"beta_2\", mu=0, sigma=10)           #定义beta_2\n",
        "    sigma = pm.Exponential(\"sigma\", 0.6)                   #定义sigma\n",
        "\n",
        "    stress = pm.MutableData(\"stress\",df.stress, dims=\"obs_id\")                           #自变量1:压力水平\n",
        "    romantic = pm.MutableData(\"romantic\",df.romantic, dims=\"obs_id\")                     #自变量2:恋爱情况\n",
        "    mu = pm.Deterministic(\"mu\", beta_0 + beta_1*stress + beta_2*romantic, dims=\"obs_id\") #定义mu，自变量与先验结合\n",
        "\n",
        "    likelihood = pm.Normal(\"y_est\", mu=mu, sigma=sigma, observed=df.networksize, dims=\"obs_id\") "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "collapsed": false,
        "id": "FA42503B631C43CC94D2C2E4CAC0D383",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/FA42503B631C43CC94D2C2E4CAC0D383/s4tt63ds2g.svg\">"
            ],
            "text/plain": [
              "<graphviz.graphs.Digraph at 0x7f803446a580>"
            ]
          },
          "execution_count": 4,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# 可视化模型（贝叶斯变量因果图）。 \n",
        "# 我们使用 pymc 自带的 `model_to_graphviz` 方法来可视化模型中各变量的因果关系。\n",
        "pm.model_to_graphviz(normal_model)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "collapsed": false,
        "id": "C760F8EE1B4C43AFAC4EF0BF408C0B02",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Auto-assigning NUTS sampler...\n",
            "Initializing NUTS using jitter+adapt_diag...\n",
            "Multiprocess sampling (4 chains in 4 jobs)\n",
            "NUTS: [beta_0, beta_1, beta_2, sigma]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 34 seconds.\n",
            "Sampling: [y_est]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "#===========================\n",
        "#     注意！！！以下代码可能需要运行1-2分钟左右\n",
        "#===========================\n",
        "with normal_model:\n",
        "    # MCMC 近似后验分布\n",
        "    normal_model_trace = pm.sample(\n",
        "        draws=5000,                   # 使用mcmc方法进行采样，draws为采样次数\n",
        "        tune=1000,                    # tune为调整采样策略的次数，可以决定这些结果是否要被保留\n",
        "        chains=4,                     # 链数\n",
        "        discard_tuned_samples= True,  # tune的结果将在采样结束后被丢弃\n",
        "        random_seed=84735)\n",
        "    \n",
        "    # 后验预测\n",
        "    normal_model_ppc = pm.sample_posterior_predictive(normal_model_trace, random_seed=84735)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "collapsed": false,
        "id": "7D8F478D2AF045EFBC6D0405FAF503F0",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "skip"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 定义计算 MAE 函数\n",
        "from statistics import median\n",
        "def MAE(model_ppc):\n",
        "    # 计算每个X取值下对应的后验预测模型的均值\n",
        "    pre_y = model_ppc.posterior_predictive[\"y_est\"].stack(sample=(\"chain\", \"draw\"))\n",
        "    pre_y_mean = pre_y.mean(axis=1).values\n",
        "\n",
        "    # 提取观测值Y，提取对应Y值下的后验预测模型的均值\n",
        "    MAE = pd.DataFrame({\n",
        "        \"networksize_ppc_mean\": pre_y_mean,\n",
        "        \"networksize_original\": df.networksize\n",
        "    })\n",
        "\n",
        "    # 计算预测误差\n",
        "    MAE[\"pre_error\"] = abs(MAE[\"networksize_original\"] -\\\n",
        "                            MAE[\"networksize_ppc_mean\"])\n",
        "\n",
        "    # 最后，计算预测误差的中位数\n",
        "    MAE = median(MAE.pre_error)\n",
        "    return MAE"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "collapsed": false,
        "id": "549B01CBDF984AC7AEEE4B1152286F1E",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "skip"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 定义计算超出HDI的数量\n",
        "def counter_outlier(model_ppc, hdi_prob=0.95):\n",
        "    # 将az.summary生成的结果存到hdi_multi这个变量中，该变量为数据框\n",
        "    hdi = az.summary(model_ppc, kind=\"stats\", hdi_prob=hdi_prob)\n",
        "    lower = hdi.iloc[:,2].values\n",
        "    upper = hdi.iloc[:,3].values\n",
        "\n",
        "    # 将原数据中的自我控制分数合并，便于后续进行判断\n",
        "    y_obs = model_ppc.observed_data[\"y_est\"].values\n",
        "\n",
        "    # 判断原数据中的压力分数是否在后验预测的95%可信区间内，并计数\n",
        "    hdi[\"verify\"] = (y_obs <= lower) | (y_obs >= upper)\n",
        "    hdi[\"y_obs\"] = y_obs\n",
        "    hdi_num = sum(hdi[\"verify\"])\n",
        "\n",
        "    return hdi_num"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "collapsed": false,
        "id": "81C99FF62C254493A0962A7CF00F53C8",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "正态回归模型MAE: 5.85\n",
            "正态回归模型超出95%hdi: 8.00\n"
          ]
        }
      ],
      "source": [
        "normal_model_MAE = MAE(normal_model_ppc)\n",
        "print(f\"正态回归模型MAE: {normal_model_MAE:.2f}\")\n",
        "\n",
        "normal_model_outliers = counter_outlier(normal_model_ppc)\n",
        "print(f\"正态回归模型超出95%hdi: {normal_model_outliers:.2f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8EEEADEE914A4EE2B60CDFFF29190587",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**PPC验证分析**  \n",
        "\n",
        "从下面的结果可以发现，后验预测的分布为正态分布，然而它并不能准确捕获数据分布的特征。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "collapsed": false,
        "id": "48A8C6AE0FBB49C394B1C5A97D77B9D7",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Axes: xlabel='y_est / y_est'>"
            ]
          },
          "execution_count": 9,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/48A8C6AE0FBB49C394B1C5A97D77B9D7/s4tt7z9xmn.png\">"
            ],
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "az.plot_ppc(normal_model_ppc, num_pp_samples=500)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2DF75631A30148D6ABFBB6A9C469C4C8",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## Poisson regression  \n",
        "\n",
        "> 在lec5中，我们曾介绍过泊松分布，当时使用的例子为在亚运会时，中国每天获得的金牌数量。  \n",
        "\n",
        "* **泊松分布(Poisson distribution)** 表示在给定$\\lambda$下，用来描述时间在一定的时间或空间内发生的事件个数，$\\lambda$为事件的平均发生率或事件发生的均值  \n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/s4o80dsjlv.png?imageView2/0/w/960/h/960)  \n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AA8CE27484604B4D81DD905447045DE8",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "* 在正态回归模型中，对数据的假设为 $Y_i|\\beta_0,\\beta_1,\\beta_2,\\sigma  \\stackrel{ind}{\\sim} N(\\mu_i, \\; \\sigma^2)  \\;\\; \\text{ with } \\;\\; \\mu_i = \\beta_0 + \\beta_1X_{i1} + \\beta_2X_{i2}$  \n",
        "\n",
        "    > 即因变量的均值可以表示为自变量的线性组合，且因变量符合以该均值为中心的正态分布 <br>  \n",
        "    > 其中 X1 为压力，X2 为恋爱情况(1代表恋爱，0代表没有恋爱)  \n",
        "\n",
        "\n",
        "* 类似的，我们也可以在泊松回归模型中做出这样的假设：  \n",
        "$Y_i | \\lambda_i \\stackrel{ind}{\\sim} Pois\\left(\\lambda_i \\right) \\;\\; \\text{ with } \\;\\; \\lambda_i = \\beta_0 + \\beta_1X_{i1} + \\beta_2X_{i2}$  \n",
        "\n",
        "**需要注意的是，在公式(1)中，事件发生数量或事件发生率 $\\lambda_i$ 的范围为 0~∞之间，然而 $\\beta_0 + \\beta_1X_{i1} + \\beta_2X_{i2}$ 的预测值会出现负数，这显然是不合理的**。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6D79208244AE4931A714F247FDDB74B7",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 泊松回归模型公式  \n",
        "\n",
        "为了解决这个问题，我们可以对 $\\lambda_i$ 进行对数转换  \n",
        "\n",
        "* 我们假设 $\\lambda_i$ 的范围为 0~50。  \n",
        "* 经过对数转化，$\\text{log}(\\lambda_i)$(纵坐标)的范围变为 -4~4。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "collapsed": false,
        "id": "56869E0A8B7740A4A7E76C61316ABB55",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/56869E0A8B7740A4A7E76C61316ABB55/s4tt7zgv2g.png\">"
            ],
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# 生成 lambdas 的范围为 0到50\n",
        "lambdas = np.arange(0.01,50,0.2)\n",
        "# 进行对数转换\n",
        "log_lambdas = np.log(lambdas)\n",
        "\n",
        "plt.plot(lambdas, log_lambdas)\n",
        "plt.xlabel(\"lambdas\")\n",
        "plt.ylabel(\"log_lambdas\")\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6EF6A1AF845D42C3A3C6B29EED6FD3CB",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "* 此时，$log(\\lambda_i)$既可以是正数也可以是负数  \n",
        "$$  \n",
        "Y_i | \\beta_0,\\beta_1, \\beta_2 \\stackrel{ind}{\\sim} Pois\\left(\\lambda_i \\right)  \n",
        "\\;\\;\\; \\text{ with } \\;\\;\\;  \n",
        "\\log\\left( \\lambda_i \\right) = \\beta_0 + \\beta_1 X_{i1} + \\beta_2 X_{i2}.  \n",
        "$$  \n",
        "\n",
        "* 因此，泊松回归模型可以写成：  \n",
        "$$  \n",
        "Y_i | \\beta_0,\\beta_1, \\beta_2 \\stackrel{ind}{\\sim} Pois\\left(\\lambda_i\\right) \\;\\; \\text{ with } \\;\\; \\lambda_i = e^{\\beta_0 + \\beta_1 X_{i1} + \\beta_2 X_{i2}}  \n",
        "$$  \n",
        "\n",
        "> 📍数学基础提醒：  \n",
        "> $$  \n",
        "> \\begin{equation}  \n",
        "> \\begin{split}  \n",
        "> \\log(\\lambda) & = \\beta_0 + \\beta_1 X_1 + \\cdots \\beta_p X_p \\\\  \n",
        "> \\lambda & = e^{\\beta_0 + \\beta_1 X_1 + \\cdots \\beta_p X_p}. \\\\  \n",
        "> \\end{split}  \n",
        "> \\end{equation}  \n",
        "> $$  \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "A548D078F0F64D648D6FC19DB942F1E0",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "在这里我们演示，在给定的$x$下 $\\lambda_i = e^{\\beta_0 + \\beta_1 X_{i1}}$ 和 $\\lambda_i = \\beta_0 + \\beta_1 X_{i1}$ 的差异  \n",
        "- 其中，我们假设 x 的范围为 0~20，并且 $beta_0=0$, $beta_1=0.2$。  \n",
        "- 使用 `np.exp` 可以进行逆对数转换。  \n",
        "- 可以发现，逆对数转换的预测结果是**非线性**的。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "collapsed": false,
        "id": "66BFCA09F11A48D1921DFF657FBCF01E",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/66BFCA09F11A48D1921DFF657FBCF01E/s4tt7zut4r.png\">"
            ],
            "text/plain": [
              "<Figure size 800x500 with 2 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "x = np.arange(20)\n",
        "beta0 = 0\n",
        "beta1 = 0.2\n",
        "lambda1 = beta0 + beta1*x\n",
        "lambda2 = np.exp(beta0 + beta1*x)\n",
        "\n",
        "fig, (ax1, ax2) = plt.subplots(1,2, figsize=(8,5))\n",
        "\n",
        "ax1.plot(x, lambda1)\n",
        "ax1.set_ylabel(\"$\\lambda$\")\n",
        "ax2.plot(x, lambda2)\n",
        "ax2.set_ylabel(\"$transformed(\\lambda)$\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CEE5916FB4DE42C1B605F09680B61E8E",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 公式中各参数的意义  \n",
        "\n",
        "$$  \n",
        "Y_i | \\beta_0,\\beta_1, \\beta_2 \\stackrel{ind}{\\sim} Pois\\left(\\lambda_i\\right) \\;\\; \\text{ with } \\;\\; \\lambda_i = e^{\\beta_0 + \\beta_1 X_{i1} + \\beta_2 X_{i2}}  \n",
        "$$  \n",
        "\n",
        "* $\\beta_0$  \n",
        "\n",
        "    * 当$X_1= X_2 = 0$时，$\\lambda_i = e^{\\beta_0}$，当自变量都为0时，事件的平均发生率为 $e^{\\beta_0}$  \n",
        "\n",
        "* $\\beta_1, ..., \\beta_p$  \n",
        "\n",
        "    * 假设 $X_2= 0$，$X_1= 1$ 时有 $\\lambda_1 =e^{\\beta_0 + \\beta_1}$，$X_1= 2$ 时有 $\\lambda_2 =e^{\\beta_0 + 2\\beta_1}$  \n",
        "\n",
        "    * 此时有 $\\frac{\\lambda_2}{\\lambda_1} = e^{\\beta_1}$ ， 因此可以得到 $\\beta_1 = log(\\lambda_2) - log(\\lambda_1)$ 或者 $e^{\\beta_1} = \\frac{\\lambda_{x+1}}{\\lambda_x}$  \n",
        "    \n",
        "    > 这表明：在其他自变量不变的情况下，自变量$X_1$每增加一个单位，事件的平均发生率$\\lambda_i$就会以$e^{\\beta_1}$倍增长  \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "285C8A7EAF244AC8B8D602FDE952D8B4",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 泊松回归的预设(assumptions)  \n",
        "\n",
        "1. **计数数据：** 适用于因变量是计数数据的情况，即非负整数，例如事件发生的次数、计数的数量等。  \n",
        "\n",
        "2. **独立性：** 假设每个观测值$Y_i$是相互独立的，即一个观测值不受其他观测值的影响。  \n",
        "\n",
        "4. **线性关系：** 经对数转换后的因变量 可以写作自变量的线性组合，如：$\\log\\left( \\lambda_i \\right) = \\beta_0 + \\beta_1 X_{i1} + \\beta_2 X_{i2}+ \\beta_3 X_{i3}$  \n",
        "\n",
        "3. **均值等于方差：** 对于泊松分布$Pois(\\lambda_i)$来说，其均值和方差都与事件的发生率相等(或约等)，$E(Y) = \\text{Var}(Y) = \\lambda$。  \n",
        "\n",
        "    * 需要注意的是，这并不意味着因变量本身的分布是泊松分布。在泊松回归中，我们说的是，在给定的X下，其观测到的Y需要满足泊松分布  \n",
        "\n",
        "    * 若此假设被满足，在数据中我们可以看到，因变量均值增大的同时，其变异性也在增大  \n",
        "\n",
        "    * 该假设是泊松回归的一个重要假设，当均值远大于方差时，泊松回归可能并不是最佳的模型选择，之后我们会介绍这种情况下更适合的模型选择。  \n",
        "    \n",
        "<table>  \n",
        "\t  <tr>  \n",
        "\t\t\t<td><img src=\"https://cdn.kesci.com/upload/s4nzvo4egp.png?imageView2/0/w/960/h/960\" alt=\"\" width=\"200\" height=\"150\"> </td>  \n",
        "\t\t\t<td><img src=\"https://www.bayesrulesbook.com/bookdown_files/figure-html/ch-12-pois-assumptions-1.png\" alt=\"\" width=\"400\" height=\"200\"> </td>  \n",
        "   </tr>  \n",
        "</table>  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DA2C68D0077C4443A257D50FCAE5C288",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 模型定义  \n",
        "\n",
        "**定义泊松回归模型的先验**  \n",
        "\n",
        "> 注意：定义先验时需要额外考虑对数变化对先验范围的影响。  \n",
        "\n",
        "$$  \n",
        "\\begin{equation}  \n",
        "\\begin{array}{lcrl}  \n",
        "\\text{data:} & \\hspace{.025in} & Y_i|\\beta_0,\\beta_1,\\beta_2 & \\stackrel{ind}{\\sim} Pois\\left(\\lambda_i \\right) \\;\\; \\text{ with } \\;\\;  \n",
        "\\lambda_i = e^{\\beta_0 + \\beta_1 X_{i1} + \\beta_2 X_{i2}}\\\\  \n",
        "\\text{priors:} & & \\beta_{0} & \\sim N\\left(0, 0.5^2 \\right) \\\\  \n",
        "&& \\beta_1 & \\sim N(0, 0.1^2) \\\\  \n",
        "&& \\beta_2 & \\sim N(0, 0.5^2) \\\\  \n",
        "\\end{array}  \n",
        "\\end{equation}  \n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "collapsed": false,
        "id": "114CD62252DC46EC9397390F65F4ABD9",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "with pm.Model(coords = {\"obs_id\": df.index})  as poisson_model:\n",
        "\n",
        "    beta_0 = pm.Normal(\"beta_0\", mu=0, sigma=0.5)          #定义beta_0          \n",
        "    beta_1 = pm.Normal(\"beta_1\", mu=0, sigma=0.1)          #定义beta_1\n",
        "    beta_2 = pm.Normal(\"beta_2\", mu=0, sigma=0.5)          #定义beta_2\n",
        "\n",
        "    #自变量1:压力水平\n",
        "    stress = pm.MutableData(\"stress\",df.stress, dims=\"obs_id\")                        \n",
        "    #自变量2:性别\n",
        "    romantic = pm.MutableData(\"romantic\",df.romantic, dims=\"obs_id\")            \n",
        "    # 预测 lambda，自变量与先验结合\n",
        "    lam = pm.Deterministic(\n",
        "        \"lam\", \n",
        "        ##------------------------------------------------\n",
        "        #  注意，这里我们使用 pm.math.exp 对预测进行逆对数转换\n",
        "        #------------------------------------------------\n",
        "        pm.math.exp(beta_0 + beta_1*stress + beta_2*romantic), dims=\"obs_id\")\n",
        "\n",
        "    likelihood = pm.Poisson(\"y_est\", mu=lam, observed=df.networksize, dims=\"obs_id\") "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "236D34E5B4694148A599DDDF7F46C052",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### 先验预测检验  \n",
        "\n",
        "* 前面说到使用正态回归模型会导致事件的发生率$\\lambda_i$为负数，我们可以使用`pm.sample_prior_predictive`进行先验预测检验，来查看生成的$\\lambda_i$情况"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "collapsed": false,
        "id": "EEAE47B11C9F4B9D9CCC5A21F812C741",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling: [beta_0, beta_1, beta_2]\n",
            "Sampling: [beta_0, beta_1, beta_2]\n"
          ]
        }
      ],
      "source": [
        "# 对正态回归模型进行先验预测\n",
        "normal_prior =pm.sample_prior_predictive(samples=50, \n",
        "                                          model=normal_model,\n",
        "                                          var_names=[\"mu\"],\n",
        "                                          random_seed=84735)\n",
        "# 对泊松回归模型进行先验预测\n",
        "poisson_prior =pm.sample_prior_predictive(samples=50, \n",
        "                                          model=poisson_model,\n",
        "                                          var_names=[\"lam\"],\n",
        "                                          random_seed=84735)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {
        "collapsed": false,
        "id": "D312D27718564B1798B45CD2E8F82F67",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/D312D27718564B1798B45CD2E8F82F67/s4tt81gomt.png\">"
            ],
            "text/plain": [
              "<Figure size 1500x500 with 2 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "fig, ax = plt.subplots(1,2,figsize=(15, 5))\n",
        "\n",
        "#绘制正态回归模型中，预测变量为压力情况下的先验预测线性模型\n",
        "ax[0].plot(normal_prior.constant_data[\"stress\"], \n",
        "           normal_prior.prior[\"mu\"].stack(sample=(\"chain\", \"draw\")), c=\"k\", alpha=0.4)\n",
        "\n",
        "#设置坐标轴标题\n",
        "ax[0].set_title('Normal regression', fontsize=20) \n",
        "ax[0].set_xlabel('Stress', fontsize=16) \n",
        "ax[0].set_ylabel('Networksize', fontsize=16) \n",
        "\n",
        "#绘制绘制泊松回归模型中，预测变量为压力情况下的先验预测线性模型\n",
        "ax[1].plot(poisson_prior.constant_data[\"stress\"], \n",
        "           poisson_prior.prior[\"lam\"].stack(sample=(\"chain\", \"draw\")), c=\"k\", alpha=0.4)\n",
        "\n",
        "#设置坐标轴标题、y轴刻度显示\n",
        "ax[1].set_title('Possion regression', fontsize=20) \n",
        "ax[1].set_xlabel('Stress', fontsize=16) \n",
        "ax[1].set_ylabel('Networksize', fontsize=16) \n",
        "ax[1].ticklabel_format(axis='y', style='plain')\n",
        "\n",
        "sns.despine()\n",
        "plt.show()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BF020E7C31964265A393A651EC2AC032",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### 绘制模型图"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "collapsed": false,
        "id": "79EBBC75FEA24F0EAE252E5D51AA1072",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/79EBBC75FEA24F0EAE252E5D51AA1072/s4tt812xq7.svg\">"
            ],
            "text/plain": [
              "<graphviz.graphs.Digraph at 0x7f7fbea58280>"
            ]
          },
          "execution_count": 15,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "pm.model_to_graphviz(poisson_model)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7415F704B383429CAD8E6808DD57614D",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### MCMC采样"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {
        "collapsed": false,
        "id": "E294E09581A045F68BBF9827A7424D14",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Auto-assigning NUTS sampler...\n",
            "Initializing NUTS using jitter+adapt_diag...\n",
            "Multiprocess sampling (4 chains in 4 jobs)\n",
            "NUTS: [beta_0, beta_1, beta_2]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 30 seconds.\n"
          ]
        }
      ],
      "source": [
        "#===========================\n",
        "#     注意！！！以下代码可能需要运行1-2分钟左右\n",
        "#===========================\n",
        "with poisson_model:\n",
        "    poisson_model_trace = pm.sample(draws=5000,            # 使用mcmc方法进行采样，draws为采样次数\n",
        "                      tune=1000,                    # tune为调整采样策略的次数，可以决定这些结果是否要被保留\n",
        "                      chains=4,                     # 链数\n",
        "                      discard_tuned_samples= True,  # tune的结果将在采样结束后被丢弃\n",
        "                      random_seed=84735)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {
        "collapsed": false,
        "id": "9E7E87AC53834978999D319CE4DEB05A",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/9E7E87AC53834978999D319CE4DEB05A/s4tt9okqdu.png\">"
            ],
            "text/plain": [
              "<Figure size 1500x1000 with 6 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# 诊断后验分布估计情况\n",
        "az.plot_trace(poisson_model_trace, var_names=[\"beta_0\",\"beta_1\",\"beta_2\"],\n",
        "              figsize=(15,10))\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C6FCAE9302C34D7AB1AD5CFE3099CF4F",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 后验参数解释  \n",
        "\n",
        "以下的结果显示：  \n",
        "- $\\beta_0 = 3.3$，那么 $e^{\\beta_0} = 28$， 表明在 X1 和 X2 为 0时，个体社交网络大小为 27。  \n",
        "- $\\beta_1 = -0.001$， $e^{\\beta_0} = 0.99$， 表明压力数值每增加1点，社交网络大小增加 0.99；$\\beta_2 = 0.005$， $e^{\\beta_2} = 1$，表明恋爱后社交网络大小增加 1。  \n",
        "- 然而，$\\beta_1$ 和 $\\beta_2$ 的94%HDI均包括0，说明压力分数和恋爱情况不能有效预测社交网络的大小。  \n",
        "\n",
        "我们接下来将通过后验预测来印证参数分析的结果。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {
        "collapsed": false,
        "id": "B01F6CAA551C4BBF8D040B0C5A4EC8E6",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/B01F6CAA551C4BBF8D040B0C5A4EC8E6/s4tt9o6ef3.png\">"
            ],
            "text/plain": [
              "<Figure size 2208x552 with 3 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "az.plot_posterior(poisson_model_trace, var_names=[\"beta_0\",\"beta_1\",\"beta_2\"])\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {
        "collapsed": false,
        "id": "40BAA77864AA47F59DC260BB6110A7C3",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/40BAA77864AA47F59DC260BB6110A7C3/s4tt9o4v3o.png\">"
            ],
            "text/plain": [
              "<Figure size 2208x552 with 3 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# 通过 np.exp 将 beta 参数进行转换\n",
        "az.plot_posterior(poisson_model_trace, var_names=[\"beta_0\",\"beta_1\",\"beta_2\"], transform = np.exp)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0831B296B73042F49536E1DBD553C49D",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 后验预测回归模型"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {
        "collapsed": false,
        "id": "4B354508D24840CC8A154E5E2286E6F4",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling: [y_est]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "with poisson_model:\n",
        "    poisson_model_ppc = pm.sample_posterior_predictive(poisson_model_trace, random_seed=84735) "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {
        "collapsed": false,
        "id": "3A38E0BFDDD44273A59643801D113EBB",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 分别获取romantic=0及romantic=1时，压力变量对应的索引值\n",
        "single_index = xr.DataArray(df.query(\"romantic == 0\")[\"stress\"])[\"index\"].values\n",
        "romantic_index = xr.DataArray(df.query(\"romantic == 1\")[\"stress\"])[\"index\"].values"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {
        "collapsed": false,
        "id": "AF62F1EFE8CD4100B19C6461EED06F70",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 根据上一步获得的索引值，将参数的后验预测分布分成两部分\n",
        "# 一部分为吸烟情况下，参数的后验预测分布、压力分数、自我控制分数\n",
        "# 另一部分为不吸烟情况下，参数的后验预测分布、压力分数、自我控制分数\n",
        "single_para = poisson_model_trace.sel(obs_id=single_index)\n",
        "romantic_para = poisson_model_trace.sel(obs_id=romantic_index)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "metadata": {
        "collapsed": false,
        "id": "4E05B9E84CFD40B29C9A43CA15B85F9F",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/4E05B9E84CFD40B29C9A43CA15B85F9F/s4tt9qx1ka.png\">"
            ],
            "text/plain": [
              "<Figure size 1000x800 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "fig, axes = plt.subplots(figsize=(10, 8))\n",
        "\n",
        "#绘制smoke=0情况下的后验预测线性模型\n",
        "az.plot_lm(idata=single_para,\n",
        "        y= \"y_est\",\n",
        "        x= \"stress\",\n",
        "        y_model = single_para.posterior[\"lam\"],\n",
        "        kind_model = \"hdi\",\n",
        "        textsize=15,           \n",
        "        grid=False,\n",
        "        axes=axes,\n",
        "        legend=False,\n",
        "        y_model_fill_kwargs = {\"color\":\"orange\"},         # 设置 hdi 区间颜色\n",
        "        y_kwargs =  {\"color\":\"orangered\",\"alpha\":0.3}     # 设置观测值颜色\n",
        "        )\n",
        "\n",
        "#----------------------------------------#\n",
        "#绘制smoke=1情况下的后验预测线性模型\n",
        "az.plot_lm(idata=romantic_para,\n",
        "        y= \"y_est\",\n",
        "        x= \"stress\",\n",
        "        y_model = romantic_para.posterior[\"lam\"],\n",
        "        kind_model = \"hdi\",\n",
        "        textsize=15,\n",
        "        grid=False,\n",
        "        axes=axes,\n",
        "        y_model_fill_kwargs = {\"color\":\"blue\"},         # 设置 hdi 区间颜色\n",
        "        y_kwargs =  {\"color\":\"navy\"},                   # 设置观测值颜色\n",
        "        legend=False)\n",
        "\n",
        "#设置坐标轴标题\n",
        "axes.set_xlabel('Stress', fontsize=16) \n",
        "axes.set_ylabel('Networksize', fontsize=16) \n",
        "\n",
        "#设置图例\n",
        "handles = [\n",
        "        plt.Line2D([], [], color=\"orange\",\n",
        "                    label = 'romantic = no'),\n",
        "        plt.Line2D([], [], color=\"blue\",\n",
        "                    label = 'romantic = yes')]\n",
        "\n",
        "axes.legend(handles=handles)\n",
        "#----------------------------------------#\n",
        "# 设置坐标轴字体大小、总标题\n",
        "axes.tick_params(labelsize=14)  \n",
        "fig.suptitle('Posterior predictive regression models', fontsize=20)\n",
        "sns.despine()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BC32C8926D554B5A8CD16EA996BB727F",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 后验预测分布  \n",
        "\n",
        "从后验预测的结果看，泊松回归模型比正态分布模型更接近观测数据的形态。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {
        "collapsed": false,
        "id": "24DE1AEC36E940A9A46DB27F95BB03FD",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Axes: xlabel='y_est / y_est'>"
            ]
          },
          "execution_count": 24,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/24DE1AEC36E940A9A46DB27F95BB03FD/s4tt9va4qk.png\">"
            ],
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/24DE1AEC36E940A9A46DB27F95BB03FD/s4tt9xc388.png\">"
            ],
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "az.plot_ppc(poisson_model_ppc, num_pp_samples=500)\n",
        "az.plot_ppc(normal_model_ppc, num_pp_samples=500)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FB46FE4399BC438E8FD3C7C6912B8879",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 评估后验预测"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {
        "collapsed": false,
        "id": "8AAA272A64D448E591B5B43DEF385719",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "泊松回归模型MAE: 5.94\n",
            "泊松回归模型超出95%hdi: 47.00\n"
          ]
        }
      ],
      "source": [
        "# 输出结果\n",
        "poisson_model_MAE = MAE(poisson_model_ppc)\n",
        "print(f\"泊松回归模型MAE: {poisson_model_MAE:.2f}\")\n",
        "\n",
        "# 输出结果\n",
        "poisson_model_outliers = counter_outlier(poisson_model_ppc)\n",
        "print(f\"泊松回归模型超出95%hdi: {poisson_model_outliers:.2f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "D9181488B25C417997F794DF60355EBE",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 模型比较"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {
        "collapsed": false,
        "id": "8A9B30494B2A42B2A37A9C05D9B57466",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "text/html": [
              "\n",
              "    <div>\n",
              "      <progress value='20000' class='' max='20000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
              "      100.00% [20000/20000 00:00&lt;00:00]\n",
              "    </div>\n",
              "    "
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "text/html": [
              "\n",
              "    <div>\n",
              "      <progress value='20000' class='' max='20000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
              "      100.00% [20000/20000 00:00&lt;00:00]\n",
              "    </div>\n",
              "    "
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# 计算模型对数似然值以便进一步计算 ELPD. \n",
        "with normal_model:\n",
        "    pm.compute_log_likelihood(normal_model_trace)\n",
        "\n",
        "with poisson_model:\n",
        "   pm.compute_log_likelihood(poisson_model_trace)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "metadata": {
        "collapsed": false,
        "id": "155FBA7F5A004AED955DF02197340A48",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>rank</th>\n",
              "      <th>elpd_loo</th>\n",
              "      <th>p_loo</th>\n",
              "      <th>elpd_diff</th>\n",
              "      <th>weight</th>\n",
              "      <th>se</th>\n",
              "      <th>dse</th>\n",
              "      <th>warning</th>\n",
              "      <th>scale</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>model1(normal)</th>\n",
              "      <td>0</td>\n",
              "      <td>-637.888936</td>\n",
              "      <td>4.328048</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.892038</td>\n",
              "      <td>11.947684</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>False</td>\n",
              "      <td>log</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>model2(poisson)</th>\n",
              "      <td>1</td>\n",
              "      <td>-709.069244</td>\n",
              "      <td>9.163981</td>\n",
              "      <td>71.180309</td>\n",
              "      <td>0.107962</td>\n",
              "      <td>30.369841</td>\n",
              "      <td>18.635034</td>\n",
              "      <td>False</td>\n",
              "      <td>log</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "                 rank    elpd_loo     p_loo  elpd_diff    weight         se  \\\n",
              "model1(normal)      0 -637.888936  4.328048   0.000000  0.892038  11.947684   \n",
              "model2(poisson)     1 -709.069244  9.163981  71.180309  0.107962  30.369841   \n",
              "\n",
              "                       dse  warning scale  \n",
              "model1(normal)    0.000000    False   log  \n",
              "model2(poisson)  18.635034    False   log  "
            ]
          },
          "execution_count": 27,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "comparison_list = {\n",
        "    \"model1(normal)\":normal_model_trace,\n",
        "    \"model2(poisson)\":poisson_model_trace,\n",
        "}\n",
        "az.compare(comparison_list)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4E45B653D7BA4F90B05CD9CAA1A72E64",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "🤔 **为什么泊松回归模型的表现反而变差了?**  \n",
        "\n",
        "* 在泊松回归中，一个重要的假设是**均值与方差相等**，但实际上，我们的数据并不满足这个假设  \n",
        "\n",
        "* 我们可以将根据恋爱状况和压力将数据进行分组，将每个组内的均值与方差进行比较  \n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/s4ok20dfbr.png?imageView2/0/w/960/h/960)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "metadata": {
        "collapsed": false,
        "id": "99CC2DB0DCD2418089D9DCE502847DA4",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th>mean</th>\n",
              "      <th>var</th>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>romantic_recode</th>\n",
              "      <th>stress_group</th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th rowspan=\"7\" valign=\"top\">no</th>\n",
              "      <th>[20,25]</th>\n",
              "      <td>23.000000</td>\n",
              "      <td>NaN</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>[25,30]</th>\n",
              "      <td>23.750000</td>\n",
              "      <td>25.071429</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>[30,35]</th>\n",
              "      <td>27.529412</td>\n",
              "      <td>106.264706</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>[35,40]</th>\n",
              "      <td>28.370370</td>\n",
              "      <td>94.396011</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>[40,45]</th>\n",
              "      <td>27.848485</td>\n",
              "      <td>114.507576</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>[45,50]</th>\n",
              "      <td>22.631579</td>\n",
              "      <td>112.578947</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>[50,55]</th>\n",
              "      <td>34.666667</td>\n",
              "      <td>142.333333</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th rowspan=\"7\" valign=\"top\">yes</th>\n",
              "      <th>[20,25]</th>\n",
              "      <td>NaN</td>\n",
              "      <td>NaN</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>[25,30]</th>\n",
              "      <td>32.666667</td>\n",
              "      <td>133.466667</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>[30,35]</th>\n",
              "      <td>25.636364</td>\n",
              "      <td>112.454545</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>[35,40]</th>\n",
              "      <td>25.173913</td>\n",
              "      <td>32.422925</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>[40,45]</th>\n",
              "      <td>28.357143</td>\n",
              "      <td>47.631868</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>[45,50]</th>\n",
              "      <td>28.625000</td>\n",
              "      <td>40.267857</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>[50,55]</th>\n",
              "      <td>19.000000</td>\n",
              "      <td>NaN</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "                                   mean         var\n",
              "romantic_recode stress_group                       \n",
              "no              [20,25]       23.000000         NaN\n",
              "                [25,30]       23.750000   25.071429\n",
              "                [30,35]       27.529412  106.264706\n",
              "                [35,40]       28.370370   94.396011\n",
              "                [40,45]       27.848485  114.507576\n",
              "                [45,50]       22.631579  112.578947\n",
              "                [50,55]       34.666667  142.333333\n",
              "yes             [20,25]             NaN         NaN\n",
              "                [25,30]       32.666667  133.466667\n",
              "                [30,35]       25.636364  112.454545\n",
              "                [35,40]       25.173913   32.422925\n",
              "                [40,45]       28.357143   47.631868\n",
              "                [45,50]       28.625000   40.267857\n",
              "                [50,55]       19.000000         NaN"
            ]
          },
          "execution_count": 28,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "#此数据中，压力的范围为20-60，可以每隔5个取为一组，如20-25为第一组，25-30为第二组\n",
        "stress_bins = range(20,60,5)\n",
        "\n",
        "#设定分组名称\n",
        "stress_labels = [f'[{i},{i+5}]' for i in stress_bins][:-1]\n",
        "\n",
        "#使用pd.cut，根据压力分数生成新的一列，为压力组别\n",
        "df[\"stress_group\"] = pd.cut(df[\"stress\"], \n",
        "                            bins=stress_bins, \n",
        "                            labels=stress_labels)\n",
        "\n",
        "#根据恋爱情况、压力组对社交网络大小进行分组，并计算每个组内的均值与方差\n",
        "df.groupby(['romantic_recode','stress_group'])['networksize'].agg(['mean','var'])\n",
        "#(注：NaN是由于分组中缺乏数据)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4AD32E67642643C1B94A6D04FBD8233B",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## Negative Binomial model  \n",
        "\n",
        "> 我们可以看到总体来说组内的方差远大于均值，在这种情况下，我们认为因变量Y是过度分散(overdispersed)的。  \n",
        "\n",
        "* 当因变量是计数数据，且过度分散，不符合泊松回归的预设时，我们可以:  \n",
        "\n",
        "    (1) 选择非泊松回归模型  \n",
        "    \n",
        "    (2) 在泊松回归模型中增加一个过度分散参数(overdispersion parameter)  \n",
        "\n",
        "* 此处我们介绍另一种非泊松回归模型，使用**负二项回归模型**来对因变量进行分析"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1B44FC2954244811A765430884078EF5",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 负二项回归模型  \n",
        "\n",
        "当$Y$为计数变量($Y \\in \\{0,1,2,\\ldots\\}$)时，可以假设其服从负二项分布：  \n",
        "\n",
        "$Y | \\mu, r \\;  \\sim \\; \\text{NegBin}(\\mu, r)$  \n",
        "\n",
        "* 其中，$\\mu$ 为平均值，$r$为过度分散参数(reciprocal dispersion parameter)  \n",
        "\n",
        "* 在负二项分布中，其预设是均值和方差不相等：  \n",
        "\n",
        "    * 当$r$非常大时，$\\text{Var}(Y) \\approx E(Y)$，此时负二项分布近似于泊松分布  \n",
        "    * 当$r$比较小时，$\\text{Var}(Y) > E(Y)$，相较于泊松分布来说，Y是过度分散的  \n",
        "$$  \n",
        "E(Y|\\mu,r) = \\mu \\;\\; \\text{ and } \\;\\; \\text{Var}(Y|\\mu,r) = \\mu + \\frac{\\mu^2}{r}.  \n",
        "$$  \n",
        "\n",
        "![Image Name](https://www.bayesrulesbook.com/bookdown_files/figure-html/negbin-ch12-1.png)  \n",
        "\n",
        "> 可以简单理解为：负二项分布中的参数$\\mu$等同于泊松分布中的参数 $\\mu$，且参数$r$控制分布的变异性  \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "30D4E72DEB8A4A3CA84BFE9874A1BD1E",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**补充知识：二项分布 vs. 负二项分布**  \n",
        "\n",
        "\n",
        "相同点：  \n",
        "- 二项分布和负二项分布都是**离散概率分布**，用于描述重复进行的独立二元试验的结果。例如抛硬币的结果、成功或失败的次数等。  \n",
        "- 二项分布和负二项分布都基于**相同的假设**，即每次试验的结果是独立的，并且成功的概率在每次试验中保持不变。  \n",
        "\n",
        "区别：  \n",
        "- 二项分布描述的是进行固定次数的独立二元试验中成功的次数。  \n",
        "\t- 它的概率质量函数给出了在给定试验次数和成功概率下，成功次数的概率分布。  \n",
        "- 负二项分布描述的是进行独立二元试验**直到达到指定数量的成功次数所需的试验次数**。  \n",
        "\t- 它的概率质量函数给出了在给定成功次数、试验次数和成功概率下，达到指定成功次数所需的试验次数的概率分布。  \n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C534455CB0174375805F43E5E7258766",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "下列示意图展示了二项分布负二项分布的区别：  \n",
        "\n",
        "- 左图表示，如果事件出现(例如，硬币朝上)的概率为 0.3，那么在3次实验中 0~3次事件出现的概率。  \n",
        "\t- 可以想象，事件不出现或者出现1次概率远高于出现3次或者2次。  \n",
        "- 右图表示，如果事件出现(例如，硬币朝上)的概率为 0.3，那么事件出现3次条件下，无数次进行实验时，事件不出现的数量以及对应的概率。  \n",
        "\t- 可以想象，如果进行10次实验，事件可能出现3次，那么事件不出现的次数是6次，并且此时的概率大于进行3次实验，事件出现3次以及事件不出现0次。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "metadata": {
        "collapsed": false,
        "id": "C2F244890A4E4D4DB760A67C030FB5CB",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/C2F244890A4E4D4DB760A67C030FB5CB/s4tta2z805.png\">"
            ],
            "text/plain": [
              "<Figure size 1200x400 with 2 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "y = np.arange(0, 30)\n",
        "k = 3\n",
        "p = 0.3\n",
        "\n",
        "fig, ax = plt.subplots(1, 2, figsize=(12, 4), sharey=True)\n",
        "\n",
        "# Binomial\n",
        "ax[0].bar(y, st.binom.pmf(y, k, p))\n",
        "ax[0].set_xticks(np.linspace(3, 30, num=10))\n",
        "ax[0].set_title(f\"Y(Binomial) = Number of k event occur\")\n",
        "\n",
        "# Negative-Binomial\n",
        "ax[1].bar(y + k, st.nbinom.pmf(y, k, p))\n",
        "ax[1].set_xticks(np.linspace(3, 30, num=10))\n",
        "ax[1].set_title(f\"Y(Negative-Binomial) = Number of trials until k event occur\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8E9F573DB70A41BE9F41378D6BB43526",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 模型定义  \n",
        "\n",
        "$$  \n",
        "\\begin{equation}  \n",
        "\\begin{array}{lcrl}  \n",
        "\\text{data:} & \\hspace{.025in} & Y_i|\\beta_0,\\beta_1,\\beta_2,r & \\stackrel{ind}{\\sim} \\text{NegBin}\\left(\\mu_i, r \\right) \\;\\; \\text{ with } \\;\\; \\mu_i= e^{\\beta_0 + \\beta_1 X_{i1} + \\beta_2 X_{i2}}\\\\  \n",
        "\\text{priors:} & & \\beta_{0}  & \\sim N\\left(0, 0.5^2\\right) \\\\  \n",
        " & & \\beta_1  & \\sim N\\left(0, 0.1^2\\right) \\\\  \n",
        " & & \\beta_2  & \\sim N\\left(0, 0.5^2\\right) \\\\  \n",
        "& & r & \\sim \\text{Exp}(1)\\\\  \n",
        "\\end{array}  \n",
        "\\end{equation}  \n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 30,
      "metadata": {
        "collapsed": false,
        "id": "5A6C1A2196A54F0BAEEAEA373B6547BA",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "with pm.Model(coords = {\"obs_id\": df.index}) as nb_model:\n",
        "\n",
        "    beta_0 = pm.Normal(\"beta_0\", mu=0, sigma=0.5)          #定义beta_0          \n",
        "    beta_1 = pm.Normal(\"beta_1\", mu=0, sigma=0.1)           #定义beta_1\n",
        "    beta_2 = pm.Normal(\"beta_2\", mu=0, sigma=0.5)           #定义beta_2\n",
        "    alpha = pm.Exponential(\"alpha\", 1)\n",
        "\n",
        "    stress = pm.MutableData(\"stress\",df.stress, dims=\"obs_id\")                           #自变量1:压力水平\n",
        "    romantic = pm.MutableData(\"romantic\",df.romantic, dims=\"obs_id\")                     #自变量2:性别\n",
        "    \n",
        "    lam = pm.Deterministic(\"lam\", \n",
        "                           pm.math.exp(beta_0 + beta_1*stress + beta_2*romantic), \n",
        "                           dims=\"obs_id\")    #定义mu，自变量与先验结合\n",
        "    \n",
        "    likelihood = pm.NegativeBinomial(\"y_est\", mu=lam, alpha=alpha, observed=df.networksize, dims=\"obs_id\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3159FC4C2C8E445697B2697BA9081315",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### 绘制模型图"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 31,
      "metadata": {
        "collapsed": false,
        "id": "3C23BD7155C6478D993F95CAD704A5F8",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/3C23BD7155C6478D993F95CAD704A5F8/s4tta2d4jk.svg\">"
            ],
            "text/plain": [
              "<graphviz.graphs.Digraph at 0x7f7fbcb55250>"
            ]
          },
          "execution_count": 31,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "pm.model_to_graphviz(nb_model)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "07F0E5240DFD45BF918124958ADD1E9A",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### MCMC采样"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "metadata": {
        "collapsed": false,
        "id": "62E752EB4A1D496B84A7B9147FCB514B",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Auto-assigning NUTS sampler...\n",
            "Initializing NUTS using jitter+adapt_diag...\n",
            "Multiprocess sampling (4 chains in 4 jobs)\n",
            "NUTS: [beta_0, beta_1, beta_2, alpha]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 47 seconds.\n"
          ]
        }
      ],
      "source": [
        "#===========================\n",
        "#     注意！！！以下代码可能需要运行1-2分钟左右\n",
        "#===========================\n",
        "with nb_model:\n",
        "    nb_model_trace = pm.sample(draws=5000,            # 使用mcmc方法进行采样，draws为采样次数\n",
        "                      tune=1000,                    # tune为调整采样策略的次数，可以决定这些结果是否要被保留\n",
        "                      chains=4,                     # 链数\n",
        "                      discard_tuned_samples= True,  # tune的结果将在采样结束后被丢弃\n",
        "                      random_seed=84735)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 33,
      "metadata": {
        "collapsed": false,
        "id": "96C205335BBE4BF2BB94B7D95ECEC1A1",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/96C205335BBE4BF2BB94B7D95ECEC1A1/s4ttcb4rxo.png\">"
            ],
            "text/plain": [
              "<Figure size 1500x1500 with 8 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "az.plot_trace(nb_model_trace, var_names=[\"beta_0\",\"beta_1\",\"beta_2\",\"alpha\"],\n",
        "              figsize=(15,15),\n",
        "              compact=False)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "metadata": {
        "collapsed": false,
        "id": "AA57D5D7A59C4DF3909FBD2A63B1F3CA",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>mean</th>\n",
              "      <th>sd</th>\n",
              "      <th>hdi_3%</th>\n",
              "      <th>hdi_97%</th>\n",
              "      <th>mcse_mean</th>\n",
              "      <th>mcse_sd</th>\n",
              "      <th>ess_bulk</th>\n",
              "      <th>ess_tail</th>\n",
              "      <th>r_hat</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>beta_0</th>\n",
              "      <td>3.069</td>\n",
              "      <td>0.166</td>\n",
              "      <td>2.765</td>\n",
              "      <td>3.388</td>\n",
              "      <td>0.002</td>\n",
              "      <td>0.001</td>\n",
              "      <td>7469.0</td>\n",
              "      <td>9088.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>beta_1</th>\n",
              "      <td>0.005</td>\n",
              "      <td>0.004</td>\n",
              "      <td>-0.002</td>\n",
              "      <td>0.013</td>\n",
              "      <td>0.000</td>\n",
              "      <td>0.000</td>\n",
              "      <td>7331.0</td>\n",
              "      <td>9121.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>beta_2</th>\n",
              "      <td>0.018</td>\n",
              "      <td>0.057</td>\n",
              "      <td>-0.088</td>\n",
              "      <td>0.124</td>\n",
              "      <td>0.000</td>\n",
              "      <td>0.000</td>\n",
              "      <td>13516.0</td>\n",
              "      <td>12164.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "         mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  \\\n",
              "beta_0  3.069  0.166   2.765    3.388      0.002    0.001    7469.0    9088.0   \n",
              "beta_1  0.005  0.004  -0.002    0.013      0.000    0.000    7331.0    9121.0   \n",
              "beta_2  0.018  0.057  -0.088    0.124      0.000    0.000   13516.0   12164.0   \n",
              "\n",
              "        r_hat  \n",
              "beta_0    1.0  \n",
              "beta_1    1.0  \n",
              "beta_2    1.0  "
            ]
          },
          "execution_count": 34,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "az.summary(nb_model_trace, var_names=[\"beta_0\",\"beta_1\",\"beta_2\"])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0BD8161E454C4E129DB2411210237192",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 后验参数解释  \n",
        "\n",
        "以下的结果显示：  \n",
        "- $\\beta_0 = 3.1$，那么 $e^{\\beta_0} = 22$， 表明在 X1 和 X2 为 0时，个体社交网络大小为 22。  \n",
        "- $\\beta_1 = 0.005$， $e^{\\beta_0} = 1$， 表明压力数值每增加1点，社交网络大小增加 1；$\\beta_2 = 0.018$， $e^{\\beta_2} = 1$，表明恋爱后社交网络大小增加 1。  \n",
        "- 然而，$\\beta_1$ 和 $\\beta_2$ 的94%HDI均包括0，说明压力分数和恋爱情况不能有效预测社交网络的大小。  \n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 35,
      "metadata": {
        "collapsed": false,
        "id": "1E75F3DAB8274AEF830B740F42524E26",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/1E75F3DAB8274AEF830B740F42524E26/s4ttcb6oek.png\">"
            ],
            "text/plain": [
              "<Figure size 2208x552 with 3 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "az.plot_posterior(nb_model_trace, var_names=[\"beta_0\",\"beta_1\",\"beta_2\"])\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 36,
      "metadata": {
        "collapsed": false,
        "id": "DE9A93162DE84C39936F320CC9092D54",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/DE9A93162DE84C39936F320CC9092D54/s4ttcbs746.png\">"
            ],
            "text/plain": [
              "<Figure size 2208x552 with 3 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# 通过 np.exp 将 beta 参数进行转换\n",
        "az.plot_posterior(nb_model_trace, var_names=[\"beta_0\",\"beta_1\",\"beta_2\"], transform = np.exp)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1081A83AF88D40C08784C0263D0DEC66",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### 后验预测分布  \n",
        "\n",
        "从可视化结果可以看出，负二项分布模型对于因变量的预测优于泊松分布。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 37,
      "metadata": {
        "collapsed": false,
        "id": "FAEC423AA2C1454A891CABF26AE8D126",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling: [y_est]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "text/plain": [
              "<Axes: xlabel='y_est / y_est'>"
            ]
          },
          "execution_count": 37,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/FAEC423AA2C1454A891CABF26AE8D126/s4ttcj1sw2.png\">"
            ],
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/FAEC423AA2C1454A891CABF26AE8D126/s4ttcj533a.png\">"
            ],
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "with nb_model:\n",
        "    nb_model_ppc = pm.sample_posterior_predictive(nb_model_trace, random_seed=84735) \n",
        "\n",
        "az.plot_ppc(nb_model_ppc, num_pp_samples=500)\n",
        "az.plot_ppc(poisson_model_ppc, num_pp_samples=500)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "92CB6BF855674E749CDAC1D781517A6E",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### 评估后验预测"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 38,
      "metadata": {
        "collapsed": false,
        "id": "EEEB87FDEEED4CB1ADE0943F6D106CB3",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "负二项回归模型MAE: 6.39\n",
            "负二项回归模型超出95%hdi: 9.00\n"
          ]
        }
      ],
      "source": [
        "# 输出结果\n",
        "nb_model_MAE = MAE(nb_model_ppc)\n",
        "print(f\"负二项回归模型MAE: {nb_model_MAE:.2f}\")\n",
        "\n",
        "# 输出结果\n",
        "nb_model_outliers = counter_outlier(nb_model_ppc)\n",
        "print(f\"负二项回归模型超出95%hdi: {nb_model_outliers:.2f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "94AD82301D9B454FAC0538A80CC78B4C",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 模型比较  \n",
        "\n",
        "* 从模型比较的结果也可以发现，负二项回归模型的 ELPD 最低，表明其结果最好。  \n",
        "* 负二项和正态模型的差异为 9.24，但这个差异和两个模型的标准误接近，说明差异并不大。  \n",
        "* MAE和后验预测区间的结果也显示：似乎正态模型的结果更优一点。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 39,
      "metadata": {
        "collapsed": false,
        "id": "C4D0542B353A436590B47B3F69243CCA",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "text/html": [
              "\n",
              "    <div>\n",
              "      <progress value='20000' class='' max='20000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
              "      100.00% [20000/20000 00:01&lt;00:00]\n",
              "    </div>\n",
              "    "
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "with nb_model:\n",
        "    pm.compute_log_likelihood(nb_model_trace)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 40,
      "metadata": {
        "collapsed": false,
        "id": "E7666F0F157F4D95972916D4FE282455",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>rank</th>\n",
              "      <th>elpd_loo</th>\n",
              "      <th>p_loo</th>\n",
              "      <th>elpd_diff</th>\n",
              "      <th>weight</th>\n",
              "      <th>se</th>\n",
              "      <th>dse</th>\n",
              "      <th>warning</th>\n",
              "      <th>scale</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>model3(negative binomial)</th>\n",
              "      <td>0</td>\n",
              "      <td>-628.642831</td>\n",
              "      <td>3.197098</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.874396</td>\n",
              "      <td>9.268887</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>False</td>\n",
              "      <td>log</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>model1(normal)</th>\n",
              "      <td>1</td>\n",
              "      <td>-637.888936</td>\n",
              "      <td>4.328048</td>\n",
              "      <td>9.246105</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>11.947684</td>\n",
              "      <td>4.457875</td>\n",
              "      <td>False</td>\n",
              "      <td>log</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>model2(poisson)</th>\n",
              "      <td>2</td>\n",
              "      <td>-709.069244</td>\n",
              "      <td>9.163981</td>\n",
              "      <td>80.426413</td>\n",
              "      <td>0.125604</td>\n",
              "      <td>30.369841</td>\n",
              "      <td>21.673498</td>\n",
              "      <td>False</td>\n",
              "      <td>log</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "                           rank    elpd_loo     p_loo  elpd_diff    weight  \\\n",
              "model3(negative binomial)     0 -628.642831  3.197098   0.000000  0.874396   \n",
              "model1(normal)                1 -637.888936  4.328048   9.246105  0.000000   \n",
              "model2(poisson)               2 -709.069244  9.163981  80.426413  0.125604   \n",
              "\n",
              "                                  se        dse  warning scale  \n",
              "model3(negative binomial)   9.268887   0.000000    False   log  \n",
              "model1(normal)             11.947684   4.457875    False   log  \n",
              "model2(poisson)            30.369841  21.673498    False   log  "
            ]
          },
          "execution_count": 40,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "comparison_list = {\n",
        "    \"model1(normal)\":normal_model_trace,\n",
        "    \"model2(poisson)\":poisson_model_trace,\n",
        "    \"model3(negative binomial)\":nb_model_trace,\n",
        "}\n",
        "az.compare(comparison_list)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 41,
      "metadata": {
        "collapsed": false,
        "id": "734D144ED265482D917539228C66F2D9",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>model_name</th>\n",
              "      <th>MAE</th>\n",
              "      <th>HDI_num</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>normal_model</td>\n",
              "      <td>5.850239</td>\n",
              "      <td>5.850239</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>possion_model</td>\n",
              "      <td>5.938275</td>\n",
              "      <td>47.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>negative_binomial_model</td>\n",
              "      <td>6.387950</td>\n",
              "      <td>9.000000</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "                model_name       MAE    HDI_num\n",
              "0             normal_model  5.850239   5.850239\n",
              "1            possion_model  5.938275  47.000000\n",
              "2  negative_binomial_model  6.387950   9.000000"
            ]
          },
          "execution_count": 41,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "summary_MAE_interval = pd.DataFrame({\n",
        "    \"model_name\":[\"normal_model\",\"possion_model\",\"negative_binomial_model\"],\n",
        "    \"MAE\":[normal_model_MAE, poisson_model_MAE, nb_model_MAE],\n",
        "    \"HDI_num\":[normal_model_MAE, poisson_model_outliers, nb_model_outliers]})\n",
        "\n",
        "summary_MAE_interval"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "92A76D6546D7400B9458550FB553C2B7",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 总结  \n",
        "\n",
        "我们现在已经学习了正态、泊松和负二项分布数据结构的贝叶斯回归模型，然而还存在其他形式的回归模型，比如伽马、逆正态、对数正态等。所有这些选项都属于广义线性模型（**generalized linear model，GLM**）的一个大类。  \n",
        "* 首先，必须注意因变量 Y 的结构。Y 是离散的还是连续的？是对称的还是倾斜的？Y 的取值范围是什么？这些问题可以帮助我们确定合适的数据结构。  \n",
        "* 其次，定义数据的期望 E(Y ∣ ...)和链接函数 g(·)。  \n",
        "\t* 对于所有 GLM，E(Y ∣ ...) 对预测因子线性组合 (X1, X2, ..., Xp) 的依赖关系表示为 $g(E(Y \\mid \\ldots))=\\beta_{0}+\\beta_{1} X_{1}+\\beta_{2} X_{2}+\\cdots+\\beta_{p} X_{p}$。其中链接函数 g(⋅) 取决于数据结构。  \n",
        "\n",
        "\n",
        "| 模型     | 正态模型      | 泊松模型      | 负二项模型    |   \n",
        "| -------- | ----------------- | ---------------- | ------------------ |  \n",
        "| 数据模型 | $Y_i\\| \\beta_0, \\beta_1, \\cdots, \\beta_p, \\sigma \\sim N(\\mu_i, \\sigma^2)$ | $Y_i\\| \\beta_0, \\beta_1, \\cdots, \\beta_p \\sim Pois(\\lambda_i)$   |   $Y_i\\| \\beta_0, \\beta_1, \\cdots, \\beta_p, \\sigma \\sim NegBin(\\mu_i, r)$    |   \n",
        "| 连接函数 | $g(\\mu_i) = \\mu_i = \\beta_0 + \\beta_1 X_{i1} + \\beta_2 X_{i2} + \\cdots + \\beta_p X_{ip}$     | $g(\\lambda_i) := \\log(\\lambda_i) = \\beta_0 + \\beta_1 X_{i1} + \\beta_2 X_{i2} + \\cdots + \\beta_p X_{ip}$ |  $g(\\mu_i) := \\log(\\mu_i) = \\beta_0 + \\beta_1 X_{i1} + \\beta_2 X_{i2} + \\cdots + \\beta_p X_{ip}$  |   \n",
        "\n",
        "\n",
        "我们不需要对每一种 GLM 进行非常仔细的研究。Bambi 工具提供了各种分布形构建帮助用户构建 GLM，以此适用于不同的数据分布。  \n",
        "\n",
        "在后续的课程中，我们将开始使用 Bambi 来构建不同的贝叶斯模型。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "584301F327C64D3098575A3AE202DC13",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### 模型定义"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 42,
      "metadata": {
        "collapsed": false,
        "id": "C0ADA9746E174A10ACC74237C0E62CC3",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "       Formula: networksize ~ stress + romantic\n",
              "        Family: negativebinomial\n",
              "          Link: mu = log\n",
              "  Observations: 174\n",
              "        Priors: \n",
              "    target = mu\n",
              "        Common-level effects\n",
              "            Intercept ~ Normal(mu: 0.0, sigma: 15.4284)\n",
              "            stress ~ Normal(mu: 0.0, sigma: 0.3791)\n",
              "            romantic ~ Normal(mu: 0.0, sigma: 5.168)\n",
              "        \n",
              "        Auxiliary parameters\n",
              "            alpha ~ HalfCauchy(beta: 1.0)"
            ]
          },
          "execution_count": 42,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# 我们需要对先验进行定义，在定义时需要注意参数的命名，可以先只传入表达式\n",
        "# 查看在bambi中各参数对应的命名，接下来进行修改\n",
        "import bambi as bmb\n",
        "\n",
        "bambi_default = bmb.Model(\"networksize ~ stress + romantic\", \n",
        "                     df,\n",
        "                     family=\"negativebinomial\")\n",
        "bambi_default "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 43,
      "metadata": {
        "collapsed": false,
        "id": "3EEB51EB4445491E928CBB6CEB76FB04",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "       Formula: networksize ~ stress + romantic\n",
              "        Family: negativebinomial\n",
              "          Link: mu = log\n",
              "  Observations: 174\n",
              "        Priors: \n",
              "    target = mu\n",
              "        Common-level effects\n",
              "            Intercept ~ Normal(mu: 0.0, sigma: 0.5)\n",
              "            stress ~ Normal(mu: 0.0, sigma: 0.1)\n",
              "            romantic ~ Normal(mu: 0.0, sigma: 0.5)\n",
              "        \n",
              "        Auxiliary parameters\n",
              "            alpha ~ Exponential(lam: 1.0)"
            ]
          },
          "execution_count": 43,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "#定义先验并传入模型中\n",
        "intercept_prior = bmb.Prior(\"Normal\", mu=0, sigma=0.5)     #对应beta_0\n",
        "stress_prior = bmb.Prior(\"Normal\", mu=0, sigma=0.1)        #对应beta_1\n",
        "romantic_prior = bmb.Prior(\"Normal\", mu=0, sigma=0.5)      #对应beta_2\n",
        "alpha_prior = bmb.Prior(\"Exponential\", lam = 1)            #对应alpha\n",
        "\n",
        "# 将三个参数的先验定义在字典prior中，注意参数的命名\n",
        "priors = {\"Intercept\": intercept_prior, \n",
        "          \"stress\": stress_prior,\n",
        "          \"romantic\": romantic_prior,\n",
        "          \"alpha\": alpha_prior}\n",
        "\n",
        "#定义关系式，传入数据\n",
        "bambi_nb = bmb.Model(\"networksize ~ stress + romantic\", \n",
        "                     df,\n",
        "                     priors=priors, \n",
        "                     family=\"negativebinomial\")\n",
        "\n",
        "#总结对模型的设置\n",
        "bambi_nb"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 44,
      "metadata": {
        "collapsed": false,
        "id": "106378914DF44EFA86F8D1912C0C381F",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Auto-assigning NUTS sampler...\n",
            "Initializing NUTS using jitter+adapt_diag...\n",
            "Multiprocess sampling (4 chains in 4 jobs)\n",
            "NUTS: [networksize_alpha, Intercept, stress, romantic]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 14 seconds.\n"
          ]
        }
      ],
      "source": [
        "bambi_nb_trace = bambi_nb.fit(draws=5000,                   # 使用mcmc方法进行采样，draws为采样次数\n",
        "                              tune=1000,                    # tune为调整采样策略的次数，可以决定这些结果是否要被保留\n",
        "                              chains=4,\n",
        "                              random_seed=84735)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 45,
      "metadata": {
        "collapsed": false,
        "id": "B3AE36C755C64D178AB9C5158A9DC359",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/B3AE36C755C64D178AB9C5158A9DC359/s4ttdpqnf1.svg\">"
            ],
            "text/plain": [
              "<graphviz.graphs.Digraph at 0x7f7fbbd82580>"
            ]
          },
          "execution_count": 45,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "bambi_nb.graph()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E07488FD41824D279A0401448940180C",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 后验参数解释  \n",
        "\n",
        "* 与pymc中的结果相比，在具体的数值大小存在差异，但关键结果并未发生改变  \n",
        "* 关键结果如：$\\beta_1$ 和 $\\beta_2$ 的94%HDI均包括0，说明压力分数和恋爱情况不能有效预测社交网络的大小"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 46,
      "metadata": {
        "collapsed": false,
        "id": "6C7CB0E631AB4E32BAF2ADD8F638F45B",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>mean</th>\n",
              "      <th>sd</th>\n",
              "      <th>hdi_3%</th>\n",
              "      <th>hdi_97%</th>\n",
              "      <th>mcse_mean</th>\n",
              "      <th>mcse_sd</th>\n",
              "      <th>ess_bulk</th>\n",
              "      <th>ess_tail</th>\n",
              "      <th>r_hat</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>Intercept</th>\n",
              "      <td>3.424</td>\n",
              "      <td>0.170</td>\n",
              "      <td>3.105</td>\n",
              "      <td>3.742</td>\n",
              "      <td>0.001</td>\n",
              "      <td>0.001</td>\n",
              "      <td>32967.0</td>\n",
              "      <td>17535.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>stress</th>\n",
              "      <td>-0.004</td>\n",
              "      <td>0.004</td>\n",
              "      <td>-0.012</td>\n",
              "      <td>0.004</td>\n",
              "      <td>0.000</td>\n",
              "      <td>0.000</td>\n",
              "      <td>33005.0</td>\n",
              "      <td>17579.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>romantic</th>\n",
              "      <td>-0.001</td>\n",
              "      <td>0.056</td>\n",
              "      <td>-0.106</td>\n",
              "      <td>0.105</td>\n",
              "      <td>0.000</td>\n",
              "      <td>0.000</td>\n",
              "      <td>34554.0</td>\n",
              "      <td>17020.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>networksize_alpha</th>\n",
              "      <td>10.637</td>\n",
              "      <td>1.484</td>\n",
              "      <td>7.897</td>\n",
              "      <td>13.432</td>\n",
              "      <td>0.008</td>\n",
              "      <td>0.006</td>\n",
              "      <td>31153.0</td>\n",
              "      <td>16157.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "                     mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  \\\n",
              "Intercept           3.424  0.170   3.105    3.742      0.001    0.001   \n",
              "stress             -0.004  0.004  -0.012    0.004      0.000    0.000   \n",
              "romantic           -0.001  0.056  -0.106    0.105      0.000    0.000   \n",
              "networksize_alpha  10.637  1.484   7.897   13.432      0.008    0.006   \n",
              "\n",
              "                   ess_bulk  ess_tail  r_hat  \n",
              "Intercept           32967.0   17535.0    1.0  \n",
              "stress              33005.0   17579.0    1.0  \n",
              "romantic            34554.0   17020.0    1.0  \n",
              "networksize_alpha   31153.0   16157.0    1.0  "
            ]
          },
          "execution_count": 46,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "az.summary(bambi_nb_trace)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "F502CCEDF2614F39AB5B92C1C36C6B26",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 后验预测分布  \n",
        "\n",
        "* 同样地，我们通过`bambi_nb.predict`生成后验预测分布并绘制，结果与pymc是相似的  \n",
        "* 在使用`bambi_nb.predict`后，后验预测分布的结果`posterior_predictive`会自动存入`bambi_nb_trace`中"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 47,
      "metadata": {
        "collapsed": false,
        "id": "96540F1AC5CE4C4BA22E733287288722",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "bambi_nb.predict(bambi_nb_trace,kind=\"pps\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 48,
      "metadata": {
        "collapsed": false,
        "id": "6E5EF1B0E90B47939FFEC618EA8CAB8F",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Axes: xlabel='networksize / networksize'>"
            ]
          },
          "execution_count": 48,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/6E5EF1B0E90B47939FFEC618EA8CAB8F/s4ttdr81k5.png\">"
            ],
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "az.plot_ppc(bambi_nb_trace,num_pp_samples=500)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AD352B8661E345549D928BF03C990189",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 练习使用Bambi  \n",
        "\n",
        "* 我们使用bambi对负二项回归模型进行了定义与后验参数采样，在本节课中我们还介绍了正态回归模型和泊松回归模型，可以试着使用Bambi对上述模型进行定义  \n",
        "\n",
        "* 如：泊松回归模型  \n",
        "$$  \n",
        "\\begin{equation}  \n",
        "\\begin{array}{lcrl}  \n",
        "\\text{data:} & \\hspace{.025in} & Y_i|\\beta_0,\\beta_1,\\beta_2 & \\stackrel{ind}{\\sim} Pois\\left(\\lambda_i \\right) \\;\\; \\text{ with } \\;\\;  \n",
        "\\lambda_i = e^{\\beta_0 + \\beta_1 X_{i1} + \\beta_2 X_{i2}}\\\\  \n",
        "\\text{priors:} & & \\beta_{0} & \\sim N\\left(0, 0.5^2 \\right) \\\\  \n",
        "&& \\beta_1 & \\sim N(0, 0.1^2) \\\\  \n",
        "&& \\beta_2 & \\sim N(0, 0.5^2) \\\\  \n",
        "\\end{array}  \n",
        "\\end{equation}  \n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 49,
      "metadata": {
        "collapsed": false,
        "id": "7C8AA4B3A3944F498A3ABFBB3EED1404",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "import bambi as bmb\n",
        "\n",
        "##---------------------------------------------------------------------------\n",
        "#      提示：参照之前的代码与先验定义，对...中的内容进行修改\n",
        "#---------------------------------------------------------------------------\n",
        "\n",
        "#定义先验并传入模型中\n",
        "intercept_prior = bmb.Prior(...)     #对应beta_0\n",
        "stress_prior = bmb.Prior(...)        #对应beta_1\n",
        "romantic_prior = bmb.Prior(...)      #对应beta_2\n",
        "alpha_prior = bmb.Prior(...)         #对应alpha\n",
        "\n",
        "# 将三个参数的先验定义在字典prior中\n",
        "priors = {\"Intercept\": ..., \n",
        "          \"stress\": ...,\n",
        "          \"romantic\": ...,\n",
        "          \"alpha\": ...}\n",
        "\n",
        "#定义关系式，传入数据\n",
        "bambi_poisson = bmb.Model(\"...\", \n",
        "                     data = df,\n",
        "                     priors=..., \n",
        "                     family=\"...\")\n",
        "\n",
        "#总结对模型的设置\n",
        "bambi_poisson"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9789F6C6E59E4C559B838538F090159A",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### MCMC采样，后验参数估计"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "id": "EE7AFFC7A07A4B3EA1451DB563A8F08B",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "bambi_poisson_trace = (...).fit(draws=5000,                   # 使用mcmc方法进行采样，draws为采样次数\n",
        "                              tune=1000,                      # tune为调整采样策略的次数，可以决定这些结果是否要被保留\n",
        "                              chains=4,\n",
        "                              random_seed=84735)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "id": "DCA443D6357C4B6F9319BBA872560E02",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "az.summary(...)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E1ABAF8C811C4D7BBA3D499E1E2A7A50",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### 后验预测分布"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "id": "E88E1708263740059C38D1E11F9F2EB6",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "bambi_poisson.predict(...,kind=\"pps\")\n",
        "az.plot_ppc(...,num_pp_samples=500)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "A04FEED7C91942A497621E1C6633F1C3",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 🎯加分练习：  \n",
        "\n",
        "**回顾之前的研究问题**  \n",
        "\n",
        "**问题0:  哪种分布可以更好的描述社交网络变量（本节课主要内容）**   \n",
        "\n",
        "问题1：个体压力与社交网络大小的关系  \n",
        "\n",
        "问题2：恋爱情况与社交网络大小的关系  \n",
        "\n",
        "**问题3：个体压力和恋爱情况与社交网络大小的关系（本节课主要内容）**  \n",
        "\n",
        "问题4：个体压力和恋爱情况对社交网络大小的交互作用  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E7C0D10E60C149769EECA7CCBC848945",
        "jupyter": {},
        "notebookId": "6565a5db5d43a49ac72d2314",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**练习目标**  \n",
        "\n",
        "在本节课中，我们知道负二项分布和正态分布可以更好地描述压力与恋爱情况对社交网络大小的影响。  \n",
        "\n",
        "在lec11中，我们逐一探究了吸烟、压力与二者地交互对自我控制水平的影响。在本节课中，我们只演示了问题3，并没有分别探究上述各预测变量与社交网络的关系，可以结合上节课的内容，使用pymc或bambi对剩下的问题进行模型定义与分析。  \n",
        "\n",
        "* Model 1: `\"networksize ~ stress \"`  \n",
        "\n",
        "* Model 2: `\"networksize ~ romantic\"`  \n",
        "\n",
        "* Model 3: `\"networksize ~ stress + romantic\"`(本节课内容)  \n",
        "\n",
        "* Model 4:`\"networksize ~ stress + romantic + stress*romantic\"`  \n",
        "\n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/image/rkz1ehen1l.png?imageView2/0/w/960/h/960)"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.5.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}
